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Abstract 

A multiscale representation-based denoising method for spherical data contaminated with 
Poisson noise, the multiscale variance stabilizing transform on the sphere (MS-VSTS), has 
been previously proposed. This paper first extends this MS-VSTS to spherical two and one 
dimensions data (2D-1D), where the two first dimensions are longitude and latitude, and the 
third dimension is a meaningful physical index such as energy or time. We then introduce a 
novel multichannel deconvolution built upon the 2D-1D MS-VSTS, which allows us to get rid 
of both the noise and the blur introduced by the point spread function (PSF) in each energy (or 
time) band. The method is applied to simulated data from the Large Area Telescope (LAT), 
the main instrument of the Fermi Gamma-Ray Space Telescope, which detects high energy 
gamma-rays in a very wide energy range (from 20 MeV to more than 300 GeV), and whose PSF 
is strongly energy-dependent (from about 3.5° at 100 MeV to less than 0.1° at 10 GeV). 
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1. Introduction 



1.1. Literature overview 

The gamma-ray sky has been studied with unprecedented sensitivity and image capability 
thanks to t he Large Area Telesc ope (LAT), the main instrument of the Fermi Gamma-Ray Space 



Telescope (lAtwood et al.|[2009), in an energy range between 20 MeV to greater than 300 GeV. 



The detection of gamma-ray point sources is difficult for two main reasons : the Poisson noise 
and the instrument's point spread function (PSF). The Poisson noise is due to the fluctuations 
in the number of detected photons. Moreover, the effect of Poisson noise is strong because of the 
weakness of the fluxes of celestial gamma rays, especially outside the Galactic plane and far away 
from intense sources. The PSF width is strongly energy-dependent, varying from about 3.5° at 100 
MeV to less than 0.1° (68% containment) at 10 GeV. Owing to large-angle multiple scattering in 
the tracker, the PSF has broad tails, for which the 95%/68% containment ratio may be as large as 
3. 

An extensive lite ratur e exists on Poisson noise removal and the interested reader may refer to 
Schmitt et al. (20111 and Starck et al. (20101 for a thorough review. Motivated by new X-ray and 
gamma-ray data challe nges, several restoration method s have been released in astrophysics that 



are based on wavelets ( |Movit| [20091 iStarck e t al.||2009| [Fay et al.||2011[ |Schmitt et"at[|2010[ ) or 



spherical data was proposed in Schmitt et al. (2010). Starck et al. (2009) developed a denoising 



rat 

L0). 



the Bayesian machinery ([ Conrad et al.||2007[ Nlorris et al. 2010| . Wavelets have also been used 
for source detection in Fermi data QAbdo et al.||2010p, and a, first Poisson denoising algorithm for 



approach that effectively handles multichannel data acquired on a cartesian grid, and where the 
third dimension can be any physically meaningful index such as wavelength, energy, or time. While 
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traditional techniques integrate over all the third dimension in order to improve the signal-to- 
noise ratio (S/N) of the sources, their approach allows us to detect the sources while preserving 
all of their spectral information without sacrificing any the sensitivity. Nonetheless, none of these 
methods take into account the point spread function (PSF) of the instrument. Deconvolution can 
be very helpful in many situations such as source identification or flux estimation in the low energy 
bands where the resolution degrades severely. To the best of our knowledge, no general method for 
both multichannel denoising and deconvolving on the sphere has been developed in the literature. 



1.2. Contributions 

We propose a general framework for denoising and deconvolution that complies with all the 
above requirements, namely deals with 

1. multichannel data, 

2. acquired on the sphere, 

3. and contaminated by Poisson noise. 

Our approach builds upon the concept o f variance stabilization applied to the spherical wavelet 
transform coefficients (Starck et al.||2006|. This approach gives a multiscale representation of the 



Poisson data with variance-stabilized coefficients, which can be treated as if they were contaminated 
by a zero-mean white Gaussian noise. The developed algorithms are validated on simulated Fermi 
HEALPix multichannel cubes (nside = 256) with energy bands ranging from 50 MeV to 50 GeV. 



1.3. Paper organization 

The rest of the paper is organized as follows. In section [2] we present a new wavelet transform 
for multichannel data on the sphere, and we show how the variance stabilization transform can 
be introduced into the decomposition. Section [3] details the way that this transform can be used 
for Gaussian and Poisson noise removal. Section [4] describes our deconvolution algorithm on the 
sphere for both mono-channel and multichannel spherical data. In Section [5] we finally draw some 
conclusions and give possible perspectives of this work. 



1.4. Notations 

For a real discrete-time filter whose impulse response is h[i], h[i] = h[—i], i e Z is its time- 
reversed version. The discrete circular convolution product of two signals is written *, where the 
term circular represents periodic boundary conditions. The symbol 5[i] is the Kronecker delta. 

For the wavelet representation, the low-pass analysis filter is denoted h and the high-pass is 
taken as g = 5 — h throughout the paper. We denote the up-sampled version of h as h^[l] = h[l] 
if l/2 j e Z and otherwise. We define = h^' 1 * ■ ■ ■ * h n * h for j 3s 1 and = S. 

The scaling and wavelet functions used for the analysis (respectively, synthesis) are denoted tj> 
(with 0(|) = ^ k h[k]<j>(x — k),x e R and k e Z) and ip (with ^(f ) = 9\k]4>{ x ~ k), x e R and k e 
Z) (respectively, <fr and tp). 



2. Multiscale representation for multichannel spherical data with poisson noise 

2.1. Fast undecimated 2D-1D wavelet decomposition/reconstruction on the sphere 

Our goal is to analyze multichannel data acquired on a sphere with a non-isotropic two-and-one- 
dimensional (2D-1D) wavelet, where the two first dimensions are spatial (longitude and latitude) 
and the third dimension is either the time or the energy. Since the dimensions do not have the same 
physical meaning, it appears natural that the wavelet scale along the third dimension (energy or 
time) should not be connected to the spatial scale. Hence, we define the wavelet function as 

^(k e ,k v ,k t ) = V W) (fc e ,fc v )V (t) (fc t ), (1) 

where ip^^ is the spherical two dimensional (2D ) spatial wavelet an d ip^ the one-dimensional (ID) 
wavelet along the third dimension. Similarly to 
dyadic spatial scales. We build the discrete 2D- 



Starck et al] ( |2009[ ), 

D wavelet decomposition by first taking a spherical 
2D undecimated wavelet transform for each k t , followed by a ID wavelet transform for each spatial 
wavelet coefficient along the third dimension. 
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Hence, for a given multichannel data set on the sphere Y[kg, k v , k t ] and after applying first the 
2D spherical undecimated wavelet transform, we have the reconstruction formula 

Ji 

Y[k e ,k v ,k t ] = a Jx [fc e , fc t ] + ^] w h [k e ,k,p,kt], Vfc t , (2) 

where J\ is the number of spatial scales, a j x is the (spatial) approximation subband, and {wj 1 }j*=i 
are the (spatial) detail subbands. To simplify the notations in the sequel, we replace the two spatial 
indices by a single index k r , which corresponds to the pixel index. Equation Q now reads 

,71 

y[kr,k t ] = a, 7l [fc r ,fc t ] + ^ w h [k r ,kt\, Vfc t . (3) 
ii=i 

For each spatial location k r and each 2D wavelet scale j±, we then apply a ID wavelet transform 
along t on the spatial wavelet coefficients Wj^kr, •] such that 

Ji 

Wj^kr,^] = Wj lt j 2 [k r ,k t ] + 2_i w ju j 2 [k r ,k t ], \f(k r ,k t ), (4) 

i2=i 

where Ji is the number of scales along t. The approximation spatial subband aj t is processed in a 
similar way, hence yielding 

■h 

a Jl [k r ,k t ] = aj lt j 2 [k r ,k t ]+ ^] w Ju j 2 [k r ,k t ], V(/c r ,fc t ). (5) 

h=i 

Inserting Eqs. Q and ^ into pj, we obtain the 2D-1D spherical undecimated wavelet represen- 
tation of Y: 

Ji J2 Ji J2 

F[fc r ,A; t ] = aj u j 2 [k r) k t ] + ^ w ii^2 [ k r, h] + ^] J2 fc t ] + ^ X! ^'^[Av, A; t ] • ( 6 ) 

31 = 1 J2 = l 31 = 1.32 = 1 

In this expression, four kinds of coefficients can be distinguished: 

— Detail-detail coefficients (ji ^ J\ and ji < Ji): 

Wjij 2 [^, '] = (S - hm) * (/4d * aji-ltfcr, '] - ^ID * ttjitfcr, '])■ (7) 

— Approximation-detail coefficients (j\ = J\ and j% ^ Ji): 

W,7ij- 2 [At, '] = h lD^ 1) * a Ji [ fc r, '] - ^1D * a -h [ fc r, •]• (8) 

— Detail-approximation coefficients (j\ -i J\ and ji = J2): 

^i.^^r, '] = ^id 2) * a 3l _i[fc r , •] - /i^ 2) * a n [k r , •]. (9) 

— Approximation-approximation coefficients (ji = Ji and ji = J2): 

a Ji,J 2 [Av, •] = ^iD 2) *aj 1 [A:r ! -]- (10) 



2.2. Multi-scale variance stabilizing transform on the sphere (MS-VSTS) 



Schmitt et al.|p010 l proposed a multiscale variance stabilizing transform adapted for Poisson 
spherical data. This transform was dubbed the multi-scale variance stabilizing transform on the 
Sphere (MS-VSTS). The MS-VSTS is a multi-scale decomposition method designed for Poisson 
noise that is based on a variance stabilizing transform(VST). Poisson noise is indeed signal- 
dependent, which complicates its removal. The aim of a VST is to get rid of this signal-dependence 
by transforming a Poisson distribution into a Gaussian distribution of known variance. In a nutshell, 
the MS-VSTS consists of plugging a VST into a multi-scale transform-the isotropic undecimated 
wavelet transform on the sphere (IUWTS)- in order to realize (approximately) Gaussian zero-mean 
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multiscale coefficients with constant variance. The noise on coefficients can then be easily removed 
using Gaussian denoising methods such as wavelet shrinkage, as we see in the next section. 

The MS-VSTS scheme is defined recursively by inserting a (nonlinear) square-root VST into 
the IUWTS steps, that is 

TTTWT^ I a i = * a i~ 1 =^> MS-VSTS f a 3 = hj-i + ctj-! nu 
U W ib 1 di = Oj-i - a 3 (VST + IUWTS) \ d 3 = T^a^) - T 3 (a ) ' 

where Tj is the VST operator on scale j 



(12) 



with the VST con stants and that depend solely on the filter h and the scale level j. 



Zhang et al. (20081 showed that the MS-VSTS detail coefficients dj on locally homogeneous parts 
of the underlying intensity signal follow asymptotically a zero-mean normal distribution with an 
intensity-independent variance that relies only on the filter h and the current scale j. Consequently, 
for a given h, b oth the stabilized va riances and the constants b^ and can be pre-computed 
once and for all (Schmitt et al.|2010l. 



(13) 



2.3. Multichannel MS-VSTS 

We now extend the MS-VSTS machinery to the multichannel case. This amounts to plugging 
the VST into the spherical 2D-1D undecimated wavelet transform introduced in Section 2.1 This 
again gives rise to four types of coefficients that take the following forms: 

— Detail-detail coefficients {j\ < J\ and j% < J2): 

Wj u j 2 [k r , •] = (S - h 1D ) * [Tj^i^i (h < ?£~ 1) * o^-iffer, •]) - rjij 3 -i (^id~ X) * a h [k r , •] j 

— Approximation-detail coefficients (ji = J\ and ji ^ J2): 

w Jlth [k r , •] = Tj lt j a -i (h^ * ajl [k r , •]) - T Jld2 * a h [k r , •]) . (14) 

— Detail-approximation coefficients (Ji < J\ and ji = J2): 

w hJ2 [kr, •] = T n _^, h (h[^ ] * a^-iikr, •]) - T jlj2 (h[^ ] * a h [k r , ■]) . (15) 

— Approximation-approximation coefficients (ji = J\ and ji = J2): 

aji,J 2 [kr,-] = * aj^kr, ■]. (16) 

In summary, all 2D-1D wavelet coefficients {wj 1 j 2 }j 1 ^j lt j 2 ^j 2 are now stabilized, and the noise 
in all these wavelet coefficients is a zero-mean Gaussian with known variance that depends solely 
on h on the resolution levels (ji, j'2). As before, these variances can be easily tabulated. 



3. Application to multichannel denoising 

We define X to be the noiseless data and Y their observed noisy version. In the case of the 
additive zero-mean white Gaussian noise, we have Y ~ Af(X,a 2 ), and for the Poisson noise Y ~ 
V(X). The main objective behind denoising is to estimate X from Y. 



3.1. Warm-up: Gaussian noise 

We start with the simple and instructive case where the noise in Y is additive white Gaussian. 
As the spherical 2D- ID undecimated wavelet transform described in Section |2.1| is linear, the 
noise remains Gaussian in the transform domain. Therefore, the thresholding strategies that have 
been developed for wavelet Gaussian denoising can be applied to the spherical 2D-1D wavelet 
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transform. Denoting the thresholding operator as TH, the denoised 2D-1D estimate of X obtained 
by thresholding the wavelet coefficients in Eq. (JsJ) reads 

Ji J2 Ji J2 

X[k r ,k t ] = a Ju j 2 [k r ,k t ]+ 2j TH-(wj u ,j 2 [k r ,k t ])+ ^] TH (u>j lj2 [k r , k t ])+ J] J] TB.(w juh [k r ,k t ]). 

jl=l J2 = l =1.72 = 1 

(17) 

A typical choice of TH is the hard thresholding operator parametrized by the scalar threshold 
r 3* 0, i.e. 

TH(*) = {° if [ x|<T ' (18) 
[x otherwise. 

The threshold r is typically chosen to be between three and five times the noise standard deviation. 



3.2. Poisson noise 

The treatment of Poisson noise is far more complicated than its Gaussian counterpart, the main 
reason being that the noise variance is equal to its mean. This is where the MS-VSTS comes into 
play. The role of the MS-VSTS is indeed to get rid of this dependence of the variance on the mean 
by ensuring that the transformed coefficients are Gaussian with constant variance (without loss of 
generality, this variance can be assumed to be equal to 1). In other words, after the MS-VSTS, we 
are brought to a Gaussian denoising problem where standard thresholding approaches apply. 

Nevertheless, denoising is not straightforward because there is no explicit reconstruction for- 
mula available because of the form of the non-linear stabilization equations above. Formally, the 
stabilizing operators Tj 1 J2 and the convolution operators along the spatial and the third dimensions 
do not commute, even though the filter bank satisfies the exact reconstruction formula. To circum- 
vent this difficulty, we propose to solve this reconstruction problem by advocating an iterative 
reconstruction scheme. 



3.3. Iterative reconstruction 

We define W to be the transform operator associated with the 2D-1D IUWTS described in 
Section |2.1| and R to be its inverse transform. We define j\4 to be the multiresolution support, 
which is determined by the set of significant coefficients detected among WV at each scale j'2) 
and location (k r ,kt), i.e. 

■M = {(ji,j2,kr,kt)\(WY) juh [k r ,kt] is significant}. (19) 

We define M to be the orthogonal projector onto j\4, i.e. Vd 

{Md) n ^k r M = \T Y }r^ rM 6 **' (20) 

v !n,nl ' 1 \d jlth [k r ,k t ] otherwise. v ; 

Our goal is to seek a solution X that preserves the significant structures of the original data 
by reproducing exactly the same wavelet coefficients as those of the input data Y, but only on 
scales and at positions where significant coefficients have been detected. Furthermore, as Poisson 
intensity functions are positive by nature, a positivity constraint is imposed on the solution. It is 
clear that there are many solutions satisfying the positivity and multiresolution support consistency 
requirements, e.g. Y itself. Thus, our reconstruction problem based solely on these constraints is 
an ill-posed inverse problem that must be regularized. Typically, the solution in which we are 
interested must be sparse by involving the lowest budget of wavelet coefficients. Therefore, our 
reconstruction is formulated as a constrained sparsity-promoting minimization problem over the 
transform coefficients d 

minldlK subject to \ d £^,h] = (Wy) jlj2 [fc r , V(j 1( j 2> k r , k t ) e M, (2J) 



and the intensity estimate X is reconstructed as X = Hd, where d is a global minimizer of Eq. (21 ). 
We recall that || • ||i = 2^ j 2 k k t Mii,^^, &t]| is the £i-norm playing the role of regularization 
and is well-known to promote sparsity (|Donoho|2004|. This problem can b e solved efficiently using 



the hybrid steepest descent algorithm ( |Yamada||200ip ( Zhang et al.||2008 ), and requires about ten 



iterations in practice. Transposed into our context, its main steps can be summarized as follows: 
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Require: Input noisy data Y, a low-pass filter h, multiresolution support M. from the detection 
step, number of iterations A max . 
1: Initialize = M.WY. 
2: tor n = 1 to A max do 
3: S n) = MS n ~ 1 '>; 
4: dOO = WP+ (R ST^ n (d("))); 
5: Update the step (3 n = (A max - n)/(N max - 1). 
6: end for 

7: return X = Rtfl N <°**\ 



where P + is the orthogonal projector onto the positive orthant. and STg n is the entry- wise 
soft-thresholding operator with threshold f3 n , i.e. for xel, STp n (x) = max(0, 1 — /3 n /\x\)x. 

The final multichannel MS-VSTS Poisson noise removal algorithm is summarized in the follow- 
ing steps: 



Require: Input noisy data Y, a low-pass filter h, threshold level r. 
1: Multichannel spherical MS- VST: Apply the 2D- ID MS-VSTS to Y using Eqs. (pj-([l6l. 
2: Detection: Detect the coefficients that are above r (significant coefficients), anclget the mul- 
tiresolution support M.. 
3: Reconstruction: Apply the above algorithm with M. to get the denoised data X. 



3.4. Experiments 

The multichannel MS-VSTS algorithm has been applied to a simulated Fermi data set, with 
14 energy bands between 50 MeV and 1.58 GeV. Figures [T] and [2] depict the denoising results for 
two energy bands. The algorithm is able to recover most of the sources, even the faint ones, on 
each energy band. Even more importantly, the 2D- ID MS-VSTS denoising algorithm allows us to 
recover the spectral information for each spatial position, as can be seen from Figure [3] 



4. Deconvolution of spherical data with Poisson noise 

We now introduce a wavelet deconvolution approach for monochannel and multichannel data 
on the sphere with Poisson noise. The main idea underlying the method is to apply the MS-VSTS 
method described above. We first introduce the deconvolution problem and then describe how the 
MS-VSTS can be used to solve the deconvolution problem. 



4.1. Problem statement 

Many problems in signal and image processing can be cast as an inversion of a linear system 

Y = HX + e , (22) 

where X e X is the data to recover, Y e y is the degraded noisy observation, e is an additive 
noise, and H : X — » y is a bounded linear operator that is typically ill-behaved because it models 
an acquisition process that encounters loss of information. When H is the identity, it is just a 



denoising problem that can be treated with the previously described methods. Inverting Eq. (22 1 
is usually an ill-posed problem, which means that there is no unique and stable solution. 

Our objective is to remove the effect of the instrument's PSF. In our case, H is the convolution 
by a blurring kernel (i.e. PSF) operator that causes Y to lack the high frequency content of X. 
Furthermore, since the noise is Poisson, e has a variance profile HA. The problem at hand is then 
a deconvolution problem in the presence of Poisson noise. 

We therefore need to both regularize the problem and handle the Poisson statistics of the 
noise. To regularize this inversion problem and reduce the space of candidate solutions, one has to 
add some prior knowledge of the typical structure of the original data X. This prior information 
accounts for the smoothness of the solution and can range from the uniform smoothness assumption 
to more complex knowledge of the geometrical structures of X. 



Schmitt et al.\ Poisson Deconvolution on the Sphere 



7 



Simulated Fermi Poisson Intensity - Energy Band = 230 MeV - 360 MeV 




0.19 ^^^K ^^^m 2.8 Log {) 

Simulated Fermi Poisson Data - Energy Band = 320 MeV - 360 MeV 




o.o ^^^k ^^^^h s.a Lo S {) 

Denoised data - Energy band = 380 MeV - 360 GeV 




0.0 ^^^k- ^^^m H.a Log () 



Figure 1. Result of the multichannel Poisson denoising algorithm on simulated Fermi data over 
the energy band 220 MeV - 360 MeV. Top: Simulated intensity skymap. Middle: Simulated noisy 
skymap. Bottom: denoised skymap. Maps are on a logarithmic scale. 



In our LAT realistic simulations, the PSF width depends strongly on the energy from 6.9° at 
50 MeV to better than 0.1° at 10 GeV and above. Figure [4] shows the normalized profiles of the 
PSF for different energy bands. 
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Simulated Fermi Poisson Intensity — Energy Band = 589 MeV — 965 MeV 




0.11 ^^^H i^^H 3.2 Log () 



Simulated Fermi Poisson Data - Energy Band = 668 MeV - 965 MeV 




0.0 ^^^^K. 3.3 Log () 



Denoised data - Energy band = 689 MeV - 965 GeV 




0.0 ^^^K- .^^^H 3.3 Log () 



Figure 2. Result of the multichannel Poisson denoising algorithm on simulated Fermi data over 
the energy band 589 MeV - 965 MeV. Top: Simulated intensity skymap. Middle: Simulated noisy 
skymap. Bottom: denoised skymap. Maps are on a logarithmic scale. 



4.2. Monochannel Deconvolution 

We first consider the single-channel case. In the literature, several algorithms have been proposed 
to perform image deconvolution on a cartesian grid. The Richardson-Lucy algorithm is certainly the 
most famous in astrophysics. In this paper, we propose a regularized Richardson-Lucy algorithm 
to deconvolve data on the sphere data. 

The Richardson-Lucy algorithm originates from a fixed-point equation obtained by maximizing 
the Poisson likelihood with respect to X while preserving positivity. This entails a multiplicative 
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Figure 3. Spectrum of a single gamma-ray point source recovered using the multichannel MS- 
VSTS denoising algorithm. Top: Single gamma-ray source from simulated Fermi data integrated 
along the energy axis (left: simulated source; middle: noisy source; right: denoised source). Bottom: 
Spectrum of the center of the point source: intensity as a function of the energy band with 14 energy 
bands between 50 MeV and 50 GeV (black: true simulated spectrum; cyan: restored spectrum from 
our denoising algorithm. 



update rule, starting at n = and X^ ' = 1 and iterating 

x (n+i) = x {n) £ (h t (F © HX (n) )) , (23) 



where ® (respectively ©) represents the element-wise multiplication (respectively division) between 
two vectors, and H T is the transpose of H whose action on an image consists in convolving it with 
the time-reversed version of the PSF associated with H. However, it is well-known that owing to 
the lack of regularization, the Richardson-Lucy algorithm tends to amplify the noise after a few 
iterations. 

We define as the residual at iteration n 

R {n) =y_HX (n) , (24) 

where R> n > can be written as the sum of its IUWTS detail subband {dj}i<cj<:j and the last ap- 
proximation subband aj, that is, 

J 

R {n) [k r ] = aj[k r ] + dj[k r l Vfc r . (25) 
i=i 

The wavelet transform provides a means of extracting only the significant structures from the 
residual at each iteration. With increasing number of iterations, a large part of the residual becomes 
statistically insignificant. The regularized significant residual is then, for a location k r 

R^[k r ] = aj[k r ]+ £ dj[k r ], (26) 

U,k r )tM 



where A4 is the multircsolution support defined in a similar way to Eq. (19 1. The regularized 
Richardson-Lucy scheme then becomes 

X (n+1) = p ( X (n) £ r H T f( HX («) + flW)@ HI W]]] , (27) 
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Figure 4. Normalized profile of the PSF for different energy bands as a function of the angle in 
degree. Black: 50 MeV - 82 MeV. Cyan: 220 MeV - 360 MeV. Orange: 960 MeV - 1.6 GeV. Blue: 
4.2 GeV - 6.9 GeV. Green: 19 GeV - 31 GeV. 



This algorithm is similar to Murtagh et al. (19951, except that the a trous wavelet transform is 
replaced by the undecimated wavelet transform. In the next subsection, we describe how the same 
algorithm can be extended to the multi-channel case. 



4.3. Multichannel deconvolution 

As the PSF is channel-dependent, the convolution observation model is now 

Y[;k t ] =U kt X[;k t ]+e[;k t ] 

in each channel kt, where Hfc t is the (spatial) convolution operator in channel k t with known PSF. 

The same recipe as in the monochannel case applies with the notable difference that the spherical 
2D-1D MS-VSTS is used instead of its monochannel counterpart. The multichannel multiresolution 
support A4 is obtained after thresholding these coefficients. 

We now define H to be the multichannel convolutiorj^j operator, which acts on a 2D-1D mul- 
tichannel spherical data set X by applying H t to each channel X[-, k t ] independently]^] The regu- 
larized multichannel Richardson-Lucy scheme is then 



(n+i) = j 

where R^ is the regularized (significant) residual 
R in) [k r ,k t ] = a Ju j 2 [k r ,k t ] + 



X<"' <g> (H T ( (Hi'") + i? (n) ) © HI< n - 



w 



31,32 



[k r , kf] 



Ui,32,k r ,k t )eM 



(28) 



(29) 



1 Strictly speaking, this is a slight abuse of terminology since the kernel is not channel-invariant. 

2 If X were to be vectorized by stacking the channels in a long column vector, H would be a block-diagonal 
matrix whose blocks are the circulant matrices H& t . 
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4.4. Experiments 

The algorithm was applied to the seven energy bands (50 MeV-1.58 GeV) of our simulated Fermi 
data set. Figures [5] to [8] display the deconvolution results for four energy bands. Figure [9] shows the 
performance of the multichannel MS-VSTS deconvolution algorithm for a single point source. The 
deconvolution not only effectively removes the blur and recovers sharply localized point sources, 
but also allows us to restore all of the spectral information. To get a better visual impression of 
the performance of the deconvolution algorithm. Figure [10] depicts the result of the algorithm on 
a single HEALPix face covering the Galactic plane. We find that the deconvolution is remarkably 
effective : our MS-VSTS multichannel deconvolution algorithm manages to remove a large part of 
the blur introduced by the PSF. 

Software 

The software related to this paper, MRS/Poisson, and its full documentation will be included 
in the next version of ISAP (Interactive Sparse astronomical data Analysis Packages) via the ISAP 
web site 13 

5. Conclusion 

This paper extends the MS-VSTS framework to deal with monochannel deconvolution, multi- 
channel denoising and multichannel deconvolution. Unlike the monochannel MS-VSTS, the multi- 
channel MS-VSTS fully exploits the information in the 2D-1D data set and allows us to recover 
the spectral information on the sources. As the PSF strongly depends on the energy, it is very 
important to have a multichannel method for deconvolution. Multichannel deconvolution using 
MS-VSTS removes a large part of the PSF blur and significantly improves the sharpness of the 
spatial localization of point sources. 
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Schmitt et al.: Poisson Deconvolution on the Sphere 



Simulated Fermi Poisson Intensity — Energy band. = 82 MeV — 134 MeV 
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Simulated Fermi Poisson Data — Energy band = 82 MeV — 134 MeV 
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Figure 5. Result of the multichannel deconvolution algorithm for different energy bands. Top: 
Simulated (blurred) intensity skymap. Middle: Blurred and noisy skymap. Bottom: Deconvolved 
skymap. Energy band : 82 MeV - 134 MeV. Maps are on a logarithmic scale. 
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Simulated Fermi Poisson Intensity - Energy Band = 230 MeV - 360 MeV 
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Simulated Fermi Poisson Data - Energy Band = 320 MeV - 360 MeV 
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Figure 6. Result of the multichannel deconvolution algorithm for different energy bands. Top: 
Simulated (blurred) intensity skymap. Middle: Blurred and noisy skymap. Bottom: Deconvolved 
skymap. Energy band : 220 MeV - 360 MeV. Maps are on a logarithmic scale. 



14 



Schmitt et al.: Poisson Deconvolution on the Sphere 



Simulated Fermi Poisson Intensity — Energy band = 360 MeV — 589 MeV 
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Simulated Fermi Poisson Data - Energy band = 360 MeV - 539 MeV 
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Figure 7. Result of the multichannel deconvolution algorithm for different energy bands. Top: 
Simulated (blurred) intensity skymap. Middle: Blurred and noisy skymap. Bottom: Deconvolved 
skymap. Energy band : 360 MeV - 589 MeV. Maps are on a logarithmic scale. 
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Simulated Fermi Deconvolved Data — Energy band = 589 MeV — 965 MeV 
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Figure 8. Result of the multichannel deconvolution algorithm for different energy bands. Top: 
Simulated (blurred) intensity skymap. Middle: Blurred and noisy skymap. Bottom: Deconvolved 
skymap. Energy band : 589 MeV - 965 MeV. Maps are on a logarithmic scale. 
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Schmitt et al.\ Poisson Deconvolution on the Sphere 




Figure 9. Profile of a single gamma-ray point source recovered using the multichannel MS-VSTS 
deconvolution algorithm. Top: Single gamma-ray point source on simulated (blurred) Fermi data 
(energy band: 360 MeV - 589 MeV) (left: simulated blurred source; middle: blurred noisy source; 
right: deconvolved source). Bottom: Profile of the point source (cyan: simulated spectrum; black: 
restored spectrum from the deconvolved source. 
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Figure 10. View on a single HEALPix face. Result of an application of the deconvolution algorithm 
to the Galactic plane. Top Left: Simulated Fermi Poisson intensity. Top Right: Simulated Fermi 
noisy data. Bottom: Fermi data deconvolved with multichannel MS-VSTS. Energy band: 360 MeV 
- 589 MeV. Pictures are on a logarithmic scale. 



